---
output: github_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(scipen = 999)
```

## Instructions for Replication:

To replicate the tables, figures and numbers in the paper, you can either run each chunks of this readme of sequentially run the following scripts (after having installed the required packages):

- script/01_main_model.R
- script/02_main_results.R
- script/03_appendix_results.R

## Package

### Installation

Run the following script to install all the required libraries

```{r, eval = F}
if(!require("pacman", quietly = T)) install.packages("pacman")
pacman::p_load(
  # Data Manipulation & Cleaning
  tidyverse, janitor, fs, glue, scales,
  
  # Data Visualization & Plotting
  ggthemes, ggridges, ggrepel, viridis, patchwork,
  
  # Statistical Modeling & Analysis
  brms, BradleyTerry2, icr,
  
  # Utilities & Performance
  remotes, tictoc
)

# Stan for bayesian models
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
cmdstanr::install_cmdstan(overwrite = T)

# Bayesian estimation for pairwise comparison
remotes::install_github('davidissamattos/bpcs')
```


### Load packages

```{r}
pacman::p_load(
  # Data Manipulation & Cleaning
  tidyverse, janitor, fs, glue, scales,
  
  # Data Visualization & Plotting
  ggthemes, ggridges, ggrepel, viridis, patchwork,
  
  # Statistical Modeling & Analysis
  brms, BradleyTerry2, icr,
  
  # Utilities & Performance
  remotes, tictoc
)

## Colors for visualization
ger_color  <- c(
  "AfD" = "#00ADEF", 
  "CDU/CSU" = "grey60",
  "FDP" = "#ebcc34",# "#ffed00",
  "Greens" = "#46962b",
  "SPD" = "#E2001A",
  "Left" = "#8B1A1A"
)
# knitr::purl("README.Rmd", output = "main_script.R")

if(!dir.exists("model")) dir_create("model")
if(!dir.exists("plot")) dir_create("plot")

```

## Input Data

### Comparison data

The data fed into the model consist of all comparison directly made by the respondents. The dataset consists in 7 columns: 

* outcome: result of the comparison (0: mp1 < mp2, 1: mp1 > mp2, 2: mp1 == mp2)
* mp_X: ID of MP X
* mp_party_X: Party of MP X
* respondent: Respondent's ID
* respondent party: Respondent's party

```{r, eval = T}
model_dt <- read_rds("input/0_respondents_rating.rds") %>%
  glimpse
```

### MP-Level Data

To validate our measure, we use collected data on representatives from different sources: 

* pageid: MPs' ID
* name: MPs' name
* party: MPs' party
* region: MPs' region
* constituency: MPs' consituency (party list if elected on party list)
* gender: MPs' gender
* faction: MPs' faction (Sälzer 2020)
* ccs_lr: MPs' score from Comparative Candidate Survey

```{r }
mp_dt <- read_rds("input/mp_dt.rds") %>%
  glimpse()
```

## Model computation

### Main model

```{r, eval = F}
set.seed(888)
tic("Model training")
model <- bpc(
  data = model_dt, 
  player0 = 'mp_1', 
  player1 = 'mp_2', 
  result_column = 'outcome',
  model_type = "davidson-U",
  cluster = "respondent",
  solve_ties = "none",
  priors = list(prior_lambda_std = 3.0,
                prior_S_std = 3.0,
                prior_U1_std = 3.0),
  show_chain_messages = T,
  seed = 888
)
toc()

write_rds(model, "model/0_main_model.rds")

tic("Parameter extraction")
param <- get_parameters_df(model, n_eff = T, Rhat = T)
toc()
write_rds(param, "data/0_main_params.rds")

main <- param %>% 
  mutate(param_type = str_extract(Parameter, ".*?(?=\\[|$)"), 
                page_id = str_extract(Parameter, "(?<=\\[A)\\d+"), 
                respondent = str_extract(Parameter, "(?<=\\[A\\d{1,4}\\,)\\w+(?=\\])")) %>%
  glimpse()

lambdas <- main %>%
  filter(param_type == "lambda") %>%
  transmute(pageid = as.numeric(page_id), 
                   mean = -Mean, 
                   median = -Median, 
                   low = -HPD_lower, 
                   high = -HPD_higher) %>%
  left_join(mp_dt) %>%
  mutate(party = fct_reorder(party, mean)) %>%
  glimpse()

write_rds(lambdas, "data/0_mp_params.rds")

respondent_bias <- main %>%
  filter(param_type == "U1") %>% 
  transmute(respondent, 
                   pageid = as.numeric(page_id), 
                   mean = -Mean, 
                   median = -Median, 
                   low = -HPD_lower, 
                   high = -HPD_higher) %>%
  left_join(mp_dt) %>%
  mutate(party = fct_reorder(party, mean)) %>%
  glimpse

write_rds(respondent_bias, "data/0_respondent_bias.rds")
```

### Bayesian bradley terry


```{r, eval = F}
tic("Model BT")
model_bt <- bpc(
  data = filter(model_dt, outcome != 2), 
  player0 = 'mp_1', 
  player1 = 'mp_2', 
  result_column = 'outcome',
  model_type = "bt-U",
  cluster = "respondent",
  solve_ties = "none",
  priors = list(prior_lambda_std = 3.0,
                prior_S_std = 3.0,
                prior_U1_std = 3.0),
  show_chain_messages = T,
  seed = 888
)
toc()
write_rds(model_bt, "model/1_model_bt.rds")

tic("Params BT")
param_bt <- get_parameters_df(model_bt, n_eff = T, Rhat = T)
toc()
write_rds(param_bt, "data/1_param_bt.rds")
```

### Frequentist Bradley Terry Model

```{r, eval = F}
bt_data <- model_dt %>%
  select(mp_1, mp_2, outcome) %>%
  filter(outcome != 2) %>%
  mutate(outcome = ifelse(outcome == 0, "win", "loss")) %>%
  count(mp_1, mp_2, outcome)  %>%
  pivot_wider(id_cols = c(mp_1, mp_2), names_from = outcome, values_from = n)  %>%
  mutate(loss = ifelse(is.na(loss), 0, loss), 
         win = ifelse(is.na(win), 0, win)) %>%
  arrange(mp_1) %>%
  mutate(mp_1 = as.factor(mp_1)) %>%
  arrange(mp_2) %>%
  mutate(mp_2 = as.factor(mp_2)) %>%
  glimpse

tic("Frequentist Bradley-Terry")
bt_model_obj <- BTm(cbind(win, loss), mp_1, mp_2, data = bt_data)
toc()
write_rds(bt_model_obj, "model/2_frequentist_bradley_terry.rds")

bt_preds <- BTabilities(bt_model_obj) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  clean_names() %>%
  transmute(pageid = as.numeric(str_extract(rowname, "\\d+")), 
            mean_bt = -ability, 
            low_bt =  mean_bt - 1.96*s_e,
            high_bt = mean_bt + 1.96*s_e) %>%
  glimpse

write_rds(bt_preds, "data/2_frequentist_preds.rds")
```


## Results

### Load artefacts

```{r Loading artefacts for main result}
# Raw comparisons
raw_answers <- read_rds("input/0_respondents_rating.rds") %>%
  mutate(
    id1 = as.numeric(str_extract(mp_1, "\\d+")),
    id2 = as.numeric(str_extract(mp_2, "\\d+"))
  ) %>%
  glimpse


# Helpers for individual visulation of MPs
targets <- read_rds("input/targets.rds") %>%
  mutate(pageid = as.factor(pageid)) %>%
  glimpse


# Estimates for each individual MP
lambdas <- read_rds("data/0_mp_params.rds") %>%
  mutate(pageid = pageid %>%
           as.factor() %>%
           fct_reorder(mean) %>%
           fct_reorder(as.numeric(party))
  ) %>%
  glimpse

# Respondent specific biases
respondent_bias <- read_rds("data/0_respondent_bias.rds") %>%
  glimpse
```

### Computed results

### Number of comparisons

```{r}
unique_comparisons <- raw_answers %>% filter(id1 > id2)
comparisons_count_by_mp <- raw_answers %>% count(mp_1, sort = T)
sum_stats <- summary(comparisons_count_by_mp$n)
respondent_counts <- raw_answers %>% 
  summarise(n = n()/2, 
            unique_mps = length(unique(mp_1)),
            .by = c(respondent, respondent_party)) %>%
  arrange(desc(n))
sum_stats_respondent = summary(respondent_counts$n)
sum_stats_respondent_unique_mps = summary(respondent_counts$unique_mps)


print(glue("Total number of comparisons: {nrow(unique_comparisons)}"))
print(glue("{comparisons_count_by_mp$mp_1[1]} (Angela Merkel) was compared {comparisons_count_by_mp$n[1]} times"))
print(glue("On average, each MP was compared {round(sum_stats[['Mean']], 2)} times (Median: {sum_stats[['Median']]})"))
print(glue("On average, each respondent compared {round(sum_stats_respondent[['Mean']], 2)} pairs of MPs (Median: {sum_stats_respondent[['Median']]})"))
print(glue("On average, each respondent compared {round(sum_stats_respondent_unique_mps[['Mean']], 2)} unique MPs (Median: {sum_stats_respondent_unique_mps[['Median']]})"))

```

### Inter-Coder Reliability - Krippendorff's alpha

```{r}
krippalpha_dt <- raw_answers %>% 
  filter(outcome != 2) %>%
  mutate(
    id1 = as.numeric(str_extract(mp_1, "\\d+")),
    id2 = as.numeric(str_extract(mp_2, "\\d+"))
  ) %>%
  filter(id1 < id2) %>%
  mutate(id = paste(mp_1, mp_2)) %>%
  filter(n() > 2, .by = id) %>%
  select(id, respondent, outcome) %>%
  mutate(outcome = recode(outcome, "0" = -1, "1" = 1, "2" = 0))

print(glue("Krippendorff's alpha computed with {nrow(krippalpha_dt)} comparisons featuring {nrow(count(krippalpha_dt, id))} unique paris"))

krippalpha_results <- krippalpha_dt %>%
  pivot_wider(id_cols = id, names_from = respondent, values_from = outcome) %>%
  select(-id) %>%
  as.matrix() %>%
  t() %>%
  krippalpha(data = ., metric = "nominal")

print(glue("Krippendorff's alpha of {round(krippalpha_results[['alpha']], 2)}"))
```

### External dataset

```{r}
ccs_correlation <- cor(as.numeric(lambdas$ccs_lr), lambdas$mean, use = 'complete.obs')

print(glue("The faction, provided by Sältzer (2020), was available for {lambdas %>% filter(!is.na(faction)) %>% nrow()} MPs."))
print(glue("We were able to merge {lambdas %>% filter(!is.na(ccs_lr)) %>% nrow()} from the Comparative Candidate Survey based on respondents party, region and position on the list"))
print(glue("The correlation between our estimates and the Comparative Candidate Survey is {round(ccs_correlation, 3)}"))
```

### Figure 2: Ideological distribution of MPs in the 19th Bundestag by political party

```{r}
lambdas %>% 
  ggplot(aes(x = mean, y = party, fill = party)) + 
  geom_density_ridges(alpha = .7, color = "grey50", quantile_lines = T, quantiles = 2,
                      jittered_points = TRUE,
                      position = position_points_jitter(width = 0.05, height = 0),
                      point_shape = '|', point_size = 4, point_alpha = .9) +
  scale_fill_manual(values = ger_color) +
  guides(fill = "none") +
  labs(x = "Ideological Position", y = "") +
  theme_few()

ggsave(filename = "plot/2_ideological_distribution_of_parties.png", width = 8, height = 5)
```

### Figure 3: Estimated individual positions for the members of the 19th Bundestag

```{r, fig.width=12, fig.height=6}
label_dt <- lambdas %>% 
  inner_join(targets) %>%
  mutate(
    nudge = case_when(
      # Ralph Brinkhaus, Thomas L. Kemmerich, Wolfgang Kubicki, Katja Suding, Cem Özdemir, Anton Hofreiter, Claudia Roth
      pageid %in% as.numeric(c("636", "699", "399", "404", "372", "427", "406", "431")) ~ 7,
      #Philipp Amthor, Carsten Linnemann, Helge Braun, Heiko Maas, Saskia Esken, Aydan Özoğuz, Jan Korte, Bernd Riexinger, Jürgen Trittin
      pageid %in% as.numeric(c("209", "210", "11", "20", "555", "605", "508", "284", "268", "441")) ~ -7,
      TRUE ~ 7*nudge
    )
  )

lambdas %>%
  ggplot(aes(x = pageid, y = mean, color = party)) +
  geom_point(size = .5, show.legend = T) +
  scale_color_manual(values = ger_color) +
  labs(y = "Ideological Position", x = "") +
  geom_text_repel(data = label_dt, aes(label = name), nudge_y = label_dt$nudge, max.iter = 3000, size = 3, force = 10) +
  geom_errorbar(aes(ymin = low, ymax = high), size = .2, alpha = .25, show.legend = F) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.border = element_blank(),
        legend.position = "none")


ggsave(filename = "plot/3_individual_positions.png", width = 12, height = 6)
```

### Figure 4: Comparison of legislators belonging to different partisan factions and between sister parties (CDU and CSU)

```{r}
lambdas %>% 
  filter(party %in% c("Greens", "SPD", "CDU/CSU", "Left", "AfD")) %>% 
  mutate(faction = case_when(party == "Greens" & faction == "FUNDI" ~ "Greens - Fundamentalist",
                             party == "Greens" & faction == "REAL" ~ "Greens - Realist",
                             party == "SPD" & faction == "SEEH" ~ "SPD - Seeheim Circle",
                             party == "SPD" ~ "SPD - Rest of the Party",
                             party == "CDU/CSU" &  region == "Bayern" ~ "CSU",
                             party == "CDU/CSU" ~ "CDU",
                             T ~ NA_character_)) %>%
  filter(!is.na(faction)) %>%
  mutate(faction = fct_reorder(as.character(faction), mean)) %>%
  ggplot(aes(x = mean, y = faction, fill = party)) + geom_boxplot(alpha = .5, show.legend = F) + 
  scale_fill_manual(values = ger_color) +
  theme_bw() +
  labs(y = "", x = "Ideological Position")

ggsave(filename = "plot/4_comparisons_between_factions.png", width = 8, height = 5)

```

### Figure 5: Comparison of expert-based estimates with self-placement

```{r}
set.seed(888)
ccs_correlation <- cor(as.numeric(lambdas$ccs_lr), lambdas$mean, use = 'complete.obs')

lambdas %>% 
  filter(!is.na(party) & !is.na(ccs_lr)) %>% 
  mutate(mean = 10*(mean - min(mean))/(max(mean) - min(mean))) %>%
  ggplot(aes(x = mean, y = ccs_lr, color = party)) +
  geom_jitter(height = 0.25) +
  scale_color_manual(values = ger_color) +
  theme_bw() +
  labs(subtitle = glue("Correlation: {round(ccs_correlation, 3)}"), 
       x = "Ideological Position (our estimate)", 
       y = "MPs' self-estimate (CCS)", 
       color = "") +
  theme(plot.background = element_rect(fill = "white", colour = NA), 
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/5_comparison_with_ccs.png", width = 8, height = 5)

```

### Figure 6: Posterior distribution of respondent-specific bias parameters

```{r}
custom_boxplot_stats_90 <- function(x) {
  r <- quantile(x, probs = c(0.05, 0.25, 0.5, 0.75, 0.95))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

# Function to compute 95% interval statistics for the boxplot
custom_boxplot_stats_95 <- function(x) {
  r <- quantile(x, probs = c(0.025, 0.25, 0.5, 0.75, 0.975))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

respondent_bias %>% 
  mutate(user = paste("Resp ", as.character(as.numeric(as.factor(respondent)))) %>% fct_reorder(mean)) %>%
  ggplot(aes(x = mean, y = user)) +
  stat_summary(fun.data = custom_boxplot_stats_95, geom = "boxplot", fatten = 2, width = 0.6, lwd = 0.2, color = "blue", alpha = 0) +
  stat_summary(fun.data = custom_boxplot_stats_90, geom = "boxplot", fatten = 2, width = 0.6, alpha = 0.6) +
  theme_minimal() +
  coord_cartesian(xlim = c(-1.5, 1.5)) +
  labs(x = "Estimated respondent specific bias", y = "Respondents\n") +
  theme(plot.background = element_rect(fill = "white", colour = NA), 
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/6_respondent_bias.png", width = 8, height = 5)
```


## Appendix

### Load additional data for Appendix

```{r Loading artefacts for appendix}
# Timestamps for each individual comparison
comp_timestamps <- read_rds("input/comp_with_timestamps.rds") %>%
  glimpse

# MP's info
mp_dt <- read_rds("input/mp_dt.rds") %>%
  glimpse()

# Raw parameter from Bayesian Davidson model
dav_bayesian_preds <- read_rds("data/0_main_params.rds") %>%
  mutate(param_type = str_extract(Parameter, ".*?(?=\\[|$)"),
         pageid = as.numeric(str_extract(Parameter, "(?<=\\[A)\\d+"))) %>%
  filter(param_type == "lambda") %>%
  transmute(pageid, mean = -Mean, low = -HPD_lower, high = -HPD_higher, Rhat) %>%
  left_join(mp_dt, by = "pageid") %>%
  mutate(party = fct_reorder(party, mean)) %>%
  glimpse

# Raw parameter from Bayesian BT model
bt_bayesian_preds <- read_rds("data/1_param_bt.rds") %>%
  mutate(param_type = str_extract(Parameter, ".*?(?=\\[|$)"),
         pageid = as.numeric(str_extract(Parameter, "(?<=\\[A)\\d+"))) %>%
  filter(param_type == "lambda") %>%
  transmute(pageid, mean_bt = -Mean, low_bt = -HPD_lower, high_bt = -HPD_higher) %>%
  glimpse

# Raw parameter from Frequentist BT model
bt_freq_preds <- read_rds("data/2_frequentist_preds.rds") %>%
  select(pageid, mean_bt_f = mean_bt, low_bt_f = low_bt, high_bt_f = high_bt) %>%
  glimpse

merged_estimates <- dav_bayesian_preds %>%
  left_join(bt_bayesian_preds, by = join_by(pageid)) %>%
  left_join(bt_freq_preds, by = join_by(pageid)) %>%
  glimpse

# Load the simulation results
sim_results <- read_rds("input/0_simulations_results.rds") %>%
  mutate(
    n_coders = fct_reorder(factor(n_coders), n_coders),
    n_pairs = fct_reorder(factor(glue("{n_pairs} pairs")), n_pairs)
  ) %>%
  glimpse

bias_avg_dt <- sim_results %>%
  group_by(random, n_coders, n_pairs, sort_units, prob_error) %>%
  summarise(
    across(c(cor, prop_found), list(mean = mean, sd = sd), na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  mutate(random = ifelse(random, "rand", "algo")) %>%
  pivot_wider(
    id_cols = c(n_coders, n_pairs, sort_units, prob_error),
    names_from = random,
    values_from = c(cor_mean, prop_found_mean, cor_sd, prop_found_sd, n)
  ) %>%
  mutate(
    cor_diff = cor_mean_rand - cor_mean_algo,
    prop_found_diff = prop_found_mean_rand - prop_found_mean_algo,
    cor_se = sqrt(cor_sd_rand^2 / n_rand + cor_sd_algo^2 / n_algo),
    prop_found_se = sqrt(prop_found_sd_rand^2 / n_rand + prop_found_sd_algo^2 / n_algo),
    cor_low = cor_diff - 1.96 * cor_se,
    cor_high = cor_diff + 1.96 * cor_se,
    prop_found_low = prop_found_diff - 1.96 * prop_found_se,
    prop_found_high = prop_found_diff + 1.96 * prop_found_se
  ) %>%
  glimpse
```

### A - Comparison pace

```{r}
median_time_per_comp <- comp_timestamps %>%
  summarise(median_time_diff = median(time_diff, na.rm = TRUE)) %>%
  round(1)

total_time <- comp_timestamps %>%
  # Removing any interval larger than 10 minutes
  filter(time_diff < 10*60) %>%
  summarise(total_time = sum(time_diff, na.rm = TRUE) / 60, .by = user) %>%
  summarise(total_mean_time = mean(total_time)) %>%
  round(1)

print(glue("The median time required to complete a comparison was {median_time_per_comp} seconds. Completing the survey required around {total_time} minutes."))
```

### Figure A2 - Comparison over time

```{r}
comp_timestamps %>%
  mutate(
    time = (time - min(time))/60/60,
    longer = time > 24,
    time = ifelse(longer, 24, time),
    .by = user
  ) %>%
  ggplot(aes(x =  user, y = time, colour = longer, shape = longer)) +
  geom_point(show.legend = F) +
  coord_flip() +
  theme_bw() +
  scale_colour_manual(values = c("black", "red")) +
  labs(y = "Time elapsed since the first comparison (in hours)", x = "Respondent")

ggsave(filename = "plot/A2_comparaison_over_time.png", width = 8, height = 5)
```


### Results of the simulations

#### Figure B3: Raw results of the simulations

```{r}
sim_results %>%
  mutate(sort_units = ifelse(sort_units, "Sorted", "Non-sorted"),
         random = ifelse(random, "Random", "Custom algorithm")) %>%
  ggplot(aes(x = cor, prop_found, color = n_pairs)) +
  geom_point(size = .5, alpha = .4) +
  facet_grid(n_coders~random) +
  scale_color_viridis(discrete = T, direction = -1, begin = 0, end = .9) +
  theme_minimal() +
  labs(x = "Correlation with Ground Truth", y = "% of recovered MPs", color = "") +
  scale_y_continuous(labels = percent) +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/B3_sim_raw.png", width = 8, height = 5)
```

#### Figure B4: Performance of the custom algorithm against the random exploration

```{r}
bias_avg_dt %>%
  ggplot(aes(x = cor_diff, y = prop_found_diff, color = n_coders)) +
  geom_hline(yintercept = 0, linetype = 4, color = "grey70") +
  geom_vline(xintercept = 0, linetype = 4, color = "grey70") +
  geom_pointrange(aes(xmin = cor_low, xmax = cor_high), size = .03) +
  geom_errorbar(aes(ymin = prop_found_low, ymax = prop_found_high)) +
  facet_wrap(~n_pairs) +
  labs(
    x = "\nCorrelation with GT (Random - Algorithm)\n[Negative means algorithm is better correlated with ground truth.]",
    y = "% recovered units (Random - Algorithm)\n[Negative means algorithm recovers more units.]\n",
    color = "Number of simulated coders"
  ) +
  scale_color_viridis(discrete = T, direction = -1) +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/B4_sim_agg.png", width = 8, height = 5)
```


### Figure C5 - Convergence and Gelman-Rubin R-hat statistics

```{r}
merged_estimates %>%
  ggplot(aes(x = Rhat)) +
  geom_histogram(bins = 43 ) +
  theme_bw() +
  labs(x = "\nGelman-Rubin R-hat statistics", y = "Number of parameter\n")

ggsave(filename = "plot/C5_rhat.png", width = 8, height = 5)
```

### Appendix D - Davidson vs. Bradley-Terry model

```{r}
cor_dav_bt <- round(cor(merged_estimates$mean, merged_estimates$mean_bt,use = "complete.obs"), 2)
print(glue("Correlation between Davidson and Bradley-Terry model: {cor_dav_bt}"))
```

#### Figure D6: Comparison between estimates from a Davidson and a Bradley-Terry model.

```{r}
merged_estimates %>%
  ggplot(aes(x = mean,  y = mean_bt, color = party)) +
  geom_point() +
  facet_wrap(~party, scales = "free") +
  theme_minimal() +
  scale_color_manual(values = ger_color) +
  guides(color = "none") +
  labs(x = "\nDavidson model\n(including 'similar positions')", y = "\nBradley-Terry model\n(excluding 'similar positions')") +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/D6_comp_davidson_bt.png", width = 8, height = 5)
```

#### Figure D7: Width of credible intervals across models.

```{r}
merged_estimates %>%
  mutate(unc = low - high,
         unc_bt =  low_bt - high_bt) %>%
  ggplot(aes(x = unc,  y = unc_bt, color = party)) +
  geom_point() +
  facet_wrap(~party) +
  theme_minimal() +
  scale_color_manual(values = ger_color) +
  coord_equal(xlim = c(0, 10), ylim = c(0, 10)) +
  guides(color = "none") +
  geom_abline(slope = 1, linetype = 4, color = "grey60") +
  labs(x = "\nDavidson model\n(including 'similar positions')", y = "\nBradley-Terry model\n(excluding 'similar positions')") +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/D7_width_credible_intervals.png", width = 8, height = 5)
```

### Appendix D: Bayesian vs. Frequentist

#### Figure D8: Distribution of estimates across Bayesian and Frequentist frameworks

```{r}
p1 <- merged_estimates %>%
  select(pageid, mean, low, high) %>%
  ggplot(aes(x = mean)) +
  geom_histogram() +
  theme_minimal() +
  labs(y = "Davidson model\n(Bayesian)", x = "Estimates") +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

p2 <- merged_estimates %>%
  select(pageid, mean_bt_f) %>%
  ggplot(aes(x = mean_bt_f)) +
  geom_histogram() +
  theme_minimal() +
  labs(y = "Bradley-Terry model\n(Frequentist)", x = "Estimates") +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

p1/p2

ggsave(filename = "plot/D8_benchmark_freq_hist.png", width = 8, height = 5)
```

#### Figure D9: Comparison between estimates from the Davidson model and the frequentist Bradley-Terry model

```{r}
merged_estimates %>%
  ggplot(aes(x = mean,  y = mean_bt_f, color = party)) +
  geom_point() +
  facet_wrap(~party, scales = "free") +
  theme_minimal() +
  scale_color_manual(values = ger_color) +
  guides(color = "none") +
  labs(x = "\nDavidson model\n(including 'similar positions')", y = "\nBradley-Terry model\n(excluding 'similar positions')") +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        panel.background = element_rect(fill = "white", colour = NA))

ggsave(filename = "plot/D9_benchmark_freq_est.png", width = 8, height = 5)
```


### Figure E11 - Gender

```{r}
merged_estimates %>%
  drop_na(gender) %>%
  transmute(
    mean,
    party,
    gender = fct_reorder(as.character(gender), mean),
    party_sex = paste(party, gender) %>%
      fct_reorder(as.numeric(as.factor(party)))
  ) %>%
  ggplot(aes(x = party_sex, y = mean, group = party_sex, fill = party)) +
  geom_boxplot(aes(alpha = gender), show.legend = F, outlier.alpha = 0) +
  scale_fill_manual(values = ger_color) +
  theme_bw() +
  coord_flip() +
  guides(fill = F) +
  labs(y = "Ideological Position", x = "")

ggsave(filename = "plot/E11_gender.png", width = 8, height = 5)
```

### Figure 12: Direct Mandate

```{r}
merged_estimates %>%
  filter(!party %in% c("FDP")) %>%
  drop_na(constituency) %>%
  transmute(
    mean,
    party,
    direct = ifelse(constituency != "Landesliste", "Constituency", "Party list") %>%
      fct_reorder(mean),
    party_direct = paste(party, direct) %>%
      fct_reorder(as.numeric(as.factor(party)))
  ) %>%
  ggplot(aes(x = party_direct, y = mean, alpha = direct, group = party_direct, fill = party)) +
  geom_boxplot(aes(alpha = direct), show.legend = F, outlier.alpha = 0) +
  scale_fill_manual(values = ger_color) +
  scale_color_viridis(discrete = T) +
  theme_bw() +
  coord_flip() +
  guides(fill = "none") +
  labs(y = "Ideological Position", x = "", fill = "")

ggsave(filename = "plot/E12_constituency_vs_list_mps.png", width = 8, height = 5)
```

### Figure 13: East-West comparison

```{r}
ost <- c("Mecklenburg-Vorpommern", "Brandenburg", "Sachsen", "Thüringen", "Sachsen-Anhalt", "Berlin")

merged_estimates %>%
  mutate(east = ifelse(region %in% ost, "- East", "- West")) %>%
  mutate(faction = paste(party, east, glue("({n()})")), .by = c(party, east)) %>%
  filter(!is.na(faction)) %>%
  mutate(faction = fct_reorder(faction, as.numeric(party))) %>%
  ggplot(aes(x = mean, y = faction, fill = party)) +
  geom_boxplot(show.legend = F, outlier.alpha = 0) +
  scale_fill_manual(values = ger_color) +
  theme_bw() +
  labs(y = "", x = "Ideological Position")

ggsave(filename = "plot/E13_east_divide.png", width = 8, height = 5)

```

## End








